Random gap model for graphene and graphene bilayers 
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The effect of a randomly fluctuating gap, created by a random staggered potential, is studied 
in a monolayer and a bilayer of graphene. The density of states, the one-particle scattering rate 
and transport properties (diffusion coefficient and conductivity) are calculated at the neutrality 
point. All these quantities vanish at a critical value of the average staggered potential, signaling a 
continuous transition to an insulating behavior. The calculations are based on the self-consistent 
Born approximation for the one-particle scattering rate and a massless mode of the two-particle 
Green's function which is created by spontaneous symmetry breaking. Transport quantities are 
directly linked to the one-particle scattering rate. Moreover, the effect of disorder is very weak in 
the case of a monolayer but much stronger in bilayer graphene. 

PACS numbers: 81.05.Uw,71.55.Ak,72.10.Bg,73.20.Jc 



Graphene, a sheet of carbon atoms, or bilayer graphene 
are semimetals with good conducting properties [H, 11, [1] . 
In particular, the minimal conductivity at the neutral- 
ity point (NP) is very robust and almost unaffected by 
disorder or thermal fluctuations [E 0, [E 0] ■ Recent ex- 
periments with hydrogenated graphene [7( and biased bi- 
layer graphene [E [tj E3] have revealed that a staggered 
potential (SP) can be created in graphene and bilayer 
graphene which breaks the sublattice symmetry. This 
opens a gap at the Fermi energy, leading to an insulating 
behavior. With this opportunity one enters a new field, 
where one can switch between a conducting and an insu- 
lating regime of a two-dimensional material, either by a 
chemical process (e.g. oxidation or hydrogenation) or by 
applying an external electric field 11 1 . 

It is clear that the opening of a uniform gap destroys 
the metallic state immediately. This means that the 
(minimal) conductivity at the NP drops from a finite 
value of order e 2 /h directly to zero. In a realistic system, 
however, the gap may not be uniform. This means that 
locally gaps open, whereas in other regions of the sample 
there is no gap. The situation can be compared with a 
classical random network of broken and unbroken bonds. 
The conductivity of such a network is nonzero as long 
as there is a percolating cluster of unbroken bonds. In 
such a system the transition from conducting to insult- 
ing behavior is presumably a second order percolation 
transition [121 ]. 

Disorder in graphene has been the subject of a num- 
ber of recent numerical studies ID, 3]. The results can 
be summarized by the statement that chiral-symmetry 
preserving disorder provides delocalized states whereas 
a chiral-symmetry breaking potential disorder leads to 
Anderson localization, even at the NP. 

Conductivity and other transport properties in 
graphene can be evaluated by solving the Bethe-Salpeter 
equation for the average two-particle Green's function 
(Cooperon) [IE EE El EE El- Unfortunately, the 
Bethe-Salpeter equation is usually a complex matrix 



equation which is difficult to handle. Therefore, a dif- 
ferent approach will be employed here that eliminates a 
part of the complexity by focusing on continuous sym- 
metries and spontaneous symmetry breaking. This al- 
lows us to identify a (massless) diffusion mode in the 
system with a randomly fluctuating gap. Consequently, 
diffusion can only stop when the spontaneous symmetry 
breaking vanishes. It will be discussed in this paper that 
this can happen if the average SP approaches a critical 
value. Moreover, there is no drop of the conductivity but 
a continuous decay to zero, depending on the fluctuations 
of the SP. 

model: Quasiparticles in monolayer graphene (MLG) 
or bilayer graphene (BLG) are described in tight-binding 
approximation by a nearest-neighbor hopping Hamilto- 
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where the underlying structure is either a honeycomb 
lattice (MLG) or two honeycomb lattices with Bernal 
stacking (BLG). 

The sublattice symmetry of the honeycomb lattice is 
broken by a staggered potential (SP) m r which is positive 
(negative) on sublattice A (B) [2fJ, |2lJ. Such a potential 
can be the result of chemical absorption of other atoms 
(e.g. oxygen or hydrogen Q) or of an external gate volt- 
age applied to the two layers of BLG 0. Neither in 
MLG nor in BLG the potential m r and, therefore, the 
gap is uniform, because of fluctuations in the coverage 
of the MLG by additional non-carbon atoms or by the 
fact that the graphene sheets are not planar 2^, 23. 24 1. 
Deviations from the planar structure in the form of rip- 
ples cause fluctuations in the distance of the two sheets 
in BLG which results in an inhomogeneous potential m r 
along each sheet. It is assumed that the gate voltage is 
adjusted at the NP such that in average the SP is exactly 
antisymmetric: {tua) = -(tub). 

At first glance, the Hamiltonian in Eq. ([1]) is a stan- 
dard hopping Hamiltonian with random potential m, fre- 
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quently used to study the generic case of Anderson local- 
ization 25J. The dispersion, however, is special in the 
case of graphene due to the honeycomb lattice: _at low- 
energies it consists of two valleys K and K' 



18, 2 



It 

is assumed that weak disorder scatters only at small mo- 
mentum such that intervalley scattering is not relevant. 
Then each valley contributes separately to transport, and 
the contribution of the two valleys to the conductivity a 
is additive: a — ok + ■ This allows us to consider for 
the low-energy properties a Dirac-type Hamiltonian for 
each valley separately 



H = h\G\ + h 2 <72 + ma^ 



(2) 



with Pauli matrices Oj and with hj 



hj = iVj (MLG), hi = Vf-Va, h 2 = 2ViV 2 (BLG) . 

(3) 

Vj is the lattice difference operator in j (— 1,2) direc- 
tion. Within this approximation the SP m r is a ran- 
dom variable with mean value (m r ) m = rh and variance 
((m r — rh)(m r i — fh)) m = gS rtr '. The following trans- 
port calculations will be based entirely on the Hamil- 
tonian of Eq. ([2|). In particular, the average Hamilto- 
nian {H) m can be diagonalized by Fourier transforma- 
tion and is k\o~\ + k 2 o-2 + mo~3 for MLG with eigenval- 
ues Ek = ±v / m 2 + k 2 . For BGL the average Hamilto- 
nian is ( k 2 - fcfjo i + 2fcifc 2 (T2 + mas with eigenvalues 
E k = ±Vfh 2 + k 4 . 

symmetries: Transport properties are determined by 
the model properties on large scales. The latter are 
controlled by the symmetry of the Hamiltonian and of 
the corresponding one-particle Green's function G(ie) = 
(H + ie)~ 1 . In the absence of sublattice-symmetry break- 
ing (i.e. for m = 0), the Hamiltonian H = h\Oi + Y12O2 
has a continuous chiral symmetry 



H 



'He 



= H 



(4) 



with a continuous parameter a, since H anticommutes 
with (T3. The SP term ma^ breaks the continuous chiral 
symmetry. However, the behavior under transposition 
hj = —hj for MLG and hj — hj for BLG provides a 
discrete symmetry: 



H 



-G^Oj 



= H 



(5) 



where j = 1 for MLG and j — 2 for BLG. This symme- 
try is broken for the one-particle Green's function G(ie) 
by the ie term. To see whether or not the symmetry is 
recovered for e — > 0, the difference 

G(ie) + a 3 G T {ie)a 3 = G(ie) - G(-ie) = mp{E = 0) (6) 

must be evaluated, where p{E = 0) = po is the den- 
sity of states at the NP. Here the limit e — > is implic- 
itly assumed. Thus the order parameter for spontaneous 
symmetry breaking is pq. 



conductivity: The conductivity can be calculated from 
the Kubo formula. Here we focus on interband scattering 
between states of energy lu/2 and —lu/2, which is a major 
contribution to transport near the NP. The frequency- 
dependent conductivity then reads (26j 



2h 



(7) 



where \$e) is an eigenstate of H in Eq. ^ with energy 
E. In other words, the conductivity is proportional to 
a matrix element of the position operator r k (k = 1,2) 
with respect to energy eigenfunctions from the lower and 
the upper band. The matrix element ($ w /2|^|^-w/2) is 
identical with the two-particle Green's function 



J2 r l Tr 2 [Gro(-w/2 - ie)G 0r {uj/2 + ie)} 



(8) 



This indicates that transport properties are expressed by 
the two-particle Green's function G(ie)G(—ie). Each of 
the two Green's functions, G(ie) and G(— ie), can be con- 
sidered as a random variable which are correlated due to 
the common random variable m r . Their distribution is 
defined by a joint distribution function P[G(ie), G{—ie)\. 
In terms of transport theory, both Green's functions must 
be included on equal footing. This is possible by intro- 
ducing the extended Green's function 



G(ie) = 



( G(ie) 
^ G(-ie) 



H + ie 
H- ie 



In the present case one can use the symmetry transfor- 
mation of H in Eq. (|5|) to write the extended Green's 
function as 



(T \ ( H + ie 



i&j 







H 1 



<Jo 
i<7j 



This introduces an extended Hamiltonian H — 
diag(H, H T ) which is invariant under a global "rotation" 



H - e s He s =H, S = [ ,° ^ 

a <jj 



(9) 



with continuous parameters a, a', since H anticommutes 
with S. The ie term of the Green's function also breaks 
this symmetry. According to Eq. ((6|), the symmetry is 
broken spontaneously for e — > if the density of states 
Po is nonzero. Since this is a continuous symmetry, there 



is a massless mode which describes diffusion [271 ] . Sym- 
metry breaking should be studied for average quantities. 
Therefore, the average density of states must be evalu- 
ated. 

spontaneous symmetry breaking: The average one- 
particle Green's function can be calculated from the av- 
erage Hamiltonian (H) m by empl oying th e self-consistent 
Born approximation (SCBA) [lj| [ll, |2l| 



(G(ie)) m » {{H) m - 2E)- 1 = G (ir),m s ) 



(10) 
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The self-energy E is a 2 x 2 tensor due to the spinor 
structure of the quasiparticles: E = —(irjcro + m s a^)/2. 
Scattering by the random SP produces an imaginary part 
of the self-energy 77 (i.e. a one-particle scattering rate) 
and a shift m s of the average SP m (i.e., m — ► to' = 
rn + rn s ). E is determined by the self-consistent equation 



-ga 3 ((H) m - 2E) rr L <7 3 



(11) 



For simplicity, the dc limit w ~ is considered here. 
The average density of states at the NP is proportional 
to the scattering rate: po = r)/2gn. This reflects that 
scattering by the random SP creates a nonzero density 
of states at the NP. It should be noticed that the entire 
calculation of the one-particle scattering rate 77 is based 
on the average one-particle Green's function. Therefore, 
it is unrelated to the continuous symmetry of Eq. (|9j). 
On the other hand, 77 > implies spontaneous breaking 
of this symmetry. 

Eq. (jlip can also be written in terms of two equations, 
one for the one-particle scattering rate 7/ and another for 
the shift of the SP m s , as 



77 = gli], m s = -mgl/(l + gl) 



(12) 



I is a function of m and 77 and also depends on the Hamil- 
tonian. For MLG it reads with momentum cutoff A 



7 = — In 

2tt 



1 



A 2 



and for BLG 
I 



4\A/ 2 + (to + m s ) 2 



if + (to + m s ) 2 



(A ~ 00) 



(13) 



(14) 



A nonzero solution 77 requires gl = 1 in the first part of 
Eq. ([T2"]) . such that m s = —fh/2 from the second part. 
Since the integrals I are monotonically decreasing func- 
tions for large m, a real solution with gl = 1 exists only 
for I to I < m c . For both, MLG and BLG, the solutions 
read 

?7 2 = (m 2 - 77i 2 )9(m 2 - to 2 )/4 , (15) 

where the model dependence enters only through the crit- 
ical average SP m c : 



2 A 



\TeM~9 - 1 



2Ae-f/ 7r (MLG), g/2 (BLG) . (16) 



m c is much bigger for BGL (cf. Fig. 1), a result which 
indicates that the effect of disorder is much stronger in 
BLG. This is also reflected by the scattering rate at to = 
which is r\ — m c /2. 

diffusion: The average two-particle Green's function 

K~ r )(ie) = -(Tr2[G rr >(-ie)G r > r (ie)]) m 



can be evaluated from an effective field theory [27j. If 
77 > the corresponding spontaneous breaking of the 



symmetry in Eq. ([9]) creates one massless mode, which 
is related to a diffusion propagator in Fourier space: 



1 



v/g 



K q (ie) e + Dq 2 
with the diffusion coefficient 

D = 9~ E riTr 2 [G Q , r0 (-iri)G Q ^(ir,)] . (17) 



Within this approximation the matrix element of the po- 
sition operator reads 



w/2\r k \®- u /2 



2»r 



2 



1 



H K q (u/2) 



9 

gu 1 



D . 

(18) 

Using the relation between the matrix element and the 
two-particle Green's function in Eq. ©, the diffusion 
coefficient becomes D — (gn/2)($i V \r 2 .\&-i V ). Inserting 



this on the right-hand side of Eq. (|18|) gives a simple rela- 
tion between the disorder averaged matrix element of r\ 
and the corresponding matrix element without disorder: 



((^/2\r 2 k \^^ /2 ))r 



^_($^| r 2|$_^) . (19) 



This is similar to the relation of the average one-particle 
Green's function in the SCBA of Eq. ([TO]), Like in the 
latter case, the averaging process leads to a change of 
energies uj/2 — ► irj (i.e. a replacement of the frequency 
by the scattering rate). Moreover, in the relation of the 
two-particle Green's function there is an extra prefactor 
— rf 2 1 (lo j2) 2 . It is important for the transport properties, 
since the average matrix element diverges like u>~ 2 . This 



indicates that the states 1$ 



±w/2i 



are delocalized for lo 



in the presence of weak SP disorder, and localization 
increases as one goes away from the NP. Such a behavior 
was also found for bond disorder in analytic (2ri | and in 
numerical studies [l4j ]. 

After evaluating ($i^|r?|$_i^), the results for the dif- 
fusion coefficient in Eq. ()17|) and for the conductivity 
in Eqs. ([7]), (fT5)) can be summerized in the following 
expressions 



D 



ag 



V 



4-7T ?7 2 + TO 2 /4 ' 



a 

(To ~ - 



7r rj 2 + m 2 /4 h 



(20) 



where a = 1 (a = 2) for MLG (BLG). First, this result 
indicates that the physical relevant quantity is the one- 
particle scattering rate 77. The difference between MLG 
and BLG is only due to the parameter a = 1, 2 and due 
to the 77i-dependent scattering rate 77. Second, the result 
reflects a diffusive behavior as long as the scattering rate 
■q does not vanish. Eq. (Ti"5|) gives a vanishing scattering 
rate for to = m c , where the critical value m c is twice the 
scattering rate at fh = 0. Moreover, the average density 
of states at the NP is proportional to rj. Therefore, a 
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FIG. 1: upper panel: dc conductivity in units of e 2 /h for BL 
graphene (upper curve) and ML graphene (lower curve) vs. 
the average staggered potential in, calculated from Eq. (|20[) 
for g = 1 and A = 1. lower panel: critical average staggered 
potential as a function of g (variance of the staggered potential 
fluctuations) for BL graphene (upper curve) and ML graphene 
(lower curve). 



global gap opens only for m > m c . Details of the trans- 
port properties distinguish between rh — and rh ^ 0. 

fa — 0: A fluctuating SP with g not too large has no 
effect on the conductivity. This can also be understood 
from the Einstein relation a cx Dp, since the density of 
states p is proportional and the diffusion coefficient D is 
inversely proportional to the normalized scattering rate 
Such a behavior was also observed in the chiral- 
invariant case with random bond disorder which is re- 
lated to ripples [H 0, 22, [26J. The scattering rate 77 
increases with disorder strength g (cf. Fig. 1). Con- 
sequently, the density of states at the NP po = rj/2gir 
increases with g for MLG, at least for small values of 
g, whereas it is constant for BLG. On the other hand, 
the diffusion coefficient decreases with 17/77, as a result of 
increased scattering. 

m 7^ 0: The conductivity decreases with 771 and even- 
tually goes to zero at 771 = m c . This is due to two effects, 



namely the reduction of the density of states and the 
reduction of the diffusion coefficient with rh, caused by 
a fluctuating gap. Since the product of the two quan- 
tities give the conductivity in the Einstein relation, the 
conductivity also decreases. 

The only difference between MLG and BLG in our cal- 
culation is the linear (MLG) and the quadratic (BLG) 
spectrum. This has quantitative consequences for the 
conductivity: For BLG it is twice as big as for MLG at 
771 = and also decays on a larger scale for < 777, < m c , 
since the critical value is m c = g/2 for BLG, whereas it 
is m c ~ exp(— 7r/<7) for MLG. As shown in Fig. 1, the 
conductivity of MLG vanishes at much lower values of rh. 
Remarkable is the enormous difference of the scattering 
rate between the two systems at 777, = 0. As shown in Fig. 
1, 77 is practically zero for a large interval of g, whereas 
it increases linearly with g for BLG. This indicates that 
disorder has a much stronger effect in the latter. 

Our result of the random SP represents a case that 
is different from random bond disorder (i.e. with chiral 
symmetry) and random scalar potential (breaks the chi- 
ral symmetry but not the sublattice symmetry) . The for- 
mer does not localize states at the NP, whereas the latter 
has presumably always localized states, with a very large 
localization length though. In a recent paper, Zhang et 
al. suggested a Kosterlitz-Thouless (KT) transition for 
a long-range random potential The KT transition 

is a phase transition that has no spontaneous symme- 
try breaking but a single massless mode in the ordered 
phase due to (7(1) phase fluctuations. In the case of the 
random SP the situation is very different: There is spon- 
taneous symmetry breaking in the diffusive phase due to 
77 > 0. Moreover, the symmetry of the fluctuations in 
Eq. ([9]) has two components rather than one. Therefore, 
the transition to the insulating behavior due to random 
SP cannot be linked to the conventional KT transition. 

A possible experimental realization of a random gap 
was recently observed by Adam et al. [2^]. It still re- 
mains to be seen whether or not the observed transition, 
which was studied by varying the gate voltage at a fixed 
gap, can be related to a nonzero average SP. This would 
require a tuning of the gap fluctuations and measurement 
of the local density of states. 

In conclusion, the one-particle scattering rate, the den- 
sity of states, the diffusion coefficient, and the conduc- 
tivity decrease with increasing average SP 777, and vanish 
at a critical point m c . The latter is exponentially small 
for MLG but proportional to disorder strength for BLG. 
Thus the effect of disorder is much stronger in BLG. 

This project was supported by a grant from the 
Deutsche Forschungsgemeinschaft and by the Aspen Cen- 
ter for Physics. 
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